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Abstract. - We study the propagation of a harmonic perturbation of small amplitude on 
a network of coupled identical phase oscillators prepared in a state of full synchronization. 
The perturbation is externally applied to a single oscillator, and is transmitted to the other 
oscillators through coupling. Numerical results and an approximate analytical treatment, valid 
for random and ordered networks, show that the response of each oscillator is a rather well- 
defined function of its distance from the oscillator where the external perturbation is applied. 
For small distances, the system behaves as a dissipative linear medium: the perturbation 
amplitude decreases exponentially with the distance, while propagating at constant speed. We 
suggest that the pattern of interactions may be deduced from measurements of the response of 
individual oscillators to perturbations applied at different nodes of the network. 



Synchronization phenomena have recently attracted a great deal of attention [1]. The 
emergence of synchronous behaviour in large ensembles of interacting dynamical elements is 
a paradigmatic form of collective self-organization, typical of a vast class of natural systems 
ranging from complex chemical reactions to large-scale biological processes [2,3]. Many of 
these phenomena are successfully reproduced by relatively simple mathematical models. The 
most basic form of synchronized dynamics, which can be achieved in an ensemble of identical 
periodic oscillators subject to attractive coupling, is full synchronization [4]. In such state, 
the individual motions of all oscillators coincide. 

A network of fully-synchronized coupled oscillators may be viewed as an active medium 
in a highly coherent collective state. A relevant question regarding the dynamics of this 
system is how the medium responds to an external perturbation which affects its coherence. 
Due to the coupling between oscillators, the perturbation should spread through the medium. 
Moreover, if dissipative mechanisms are acting and the state of full synchronization is stable, 
the perturbation should die out as the distance from the point where it is applied becomes 
larger. The properties of this propagation phenomenon provide a dynamical characterization 
of the synchronized medium, much like the propagation of a defect in an extended system. 
This is the problem that we address in this Letter. 

Kuramoto's model of coupled phase oscillators [5] provides the basis for our model. The 
state of each oscillator is described by a phase variable fa £ [0, 2tt) which, in the absence of 
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coupling, rotates at constant frequency, fa = u>i . A network of coupled oscillators is governed 
by the equations 

N 

fa = ojj + fc Jjj sm(fa - fa), (1) 

where k is the coupling strength. Here, Jjj = 1 if oscillator i is coupled to oscillator j, 
and Jy = otherwise. The connection (or adjacency) matrix J — {Jij} is not necessarily 
symmetric. The underlying connection network is therefore a directed graph, where the links 
pointing to the node occupied by a given oscillator i start at those oscillators which enter the 
equation of motion for fa. 

If the oscillators are identical, u>i = loj for all i, j. Without generality loss, their natural 
frequencies are chosen to be u>i = 0, and the coupling strength is fixed at k = 1. To represent 
the external perturbation, one of the oscillators -say, i — 1- is also coupled to an oscillating 
force of strength e and frequency fl. Starting from eq. JJ), the equation of motion for fa can 
be written as 

JV 

fa = Jjj sin(</>j - fa) + eSu sin(fit - fa), (2) 
i=i 

where Sij is the Kronecker symbol. In the absence of external forcing (e = 0) the ensemble 
admits a fully synchronized state where fa = fa for all i, with constant fa. The stability 
of this state can be analytically proven for some specific forms of the connection matrix J . 
In particular, full synchronization is stable in the case where the number of non- vanishing 
elements in all rows of J is the same, i.e. when the number Zi — J^ of connections ending 
at each oscillator is the same, Zi = z, for all i [6]. In this case, the network is a regular directed 
graph [7] . The following study is mainly focused on this kind of network. 

A small perturbation of the fully synchronized state produces a deviation of the same order 
e as the perturbation. For e — > 0, the solution to eq. J2J) can be found by writing fa = fa +eipi, 
and expanding up to the first order in e, which yields 

N 

ipi = Jiji^Pj - ipi) + 5n exp(iQt). (3) 
i=i 

In the last term, we have replaced sin(ftt) by exp(iftt), for convenience in the mathemati- 
cal treatment. After transients have elapsed, the solution to eq. J3J has the form ^(t) = 
Ai exp(iftt). The complex amplitude A4 is obtained from the linear system 

MA = ei, (4) 

where A = (Ax, A 2 , ■ ■ -,A N ), e x = (1,0, ... ,0), and M = (z + ifi)I - J, with 1 the N x N 
identity matrix. The amplitudes are, thus, A = Mr e\. For a given connection matrix J, 
they can be found by numerically inverting M. 

Figure n shows typical results for the moduli \Aj\ and the phases (fi of the amplitudes, 
Aj = \Ai I ex.p(i(fi), in a 10 3 -oscillator random network. Each oscillator is coupled to z = 2 
neighbours, which are chosen at random from the whole ensemble (avoiding self and multiple 
connections) . The three data sets of the figure show results for the same random network and 
different perturbation frequencies SI. Each dot represents the modulus or phase of a single 
oscillator i as a function of its distance di to oscillator 1, where the external perturbation is 
applied. This distance is defined as the number of connections along the shortest (directed) 
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Fig. 1 - Moduli (a) and phases (b) of the individual amplitudes A4 in a random network [z = 2) of 
10 3 fully synchronized oscillators, as functions of the distance to the perturbed oscillator. Each dot 
corresponds to the amplitude of a single oscillator. The three data sets have been obtained for the 
same realization of the network but with different perturbation frequencies fl. Full and dotted lines 
are analytical predictions for small and large distances, respectively. 

path starting at oscillator 1 and ending at i. For this particular network, di varies between 
(for i = 1) and 15. 

The numerical results of fig. ^ have been obtained by inversion of the matrix M. in eq. J3J. 
We have verified that, as expected, these results are reproduced by numerical integration of the 
equations of motion J5J, for sufficiently small values of the perturbation amplitude (e < 10 -4 ). 

It turns out that the modulus \Aj\ has a rather well-defined dependence on the distance 
di. For small distances, it decays exponentially, and its value is essentially the same for all 
the oscillators at a given distance. As di grows, however, the values of \Ai\ become more 
dispersed as a function of di, especially for large perturbation frequencies. Moreover, the 
exponential decay breaks down, and the variation of |j4j| with the distance becomes much 
slower. The values of \Ai\ for large distances depend on the frequency f2. Numerical results 
for different realizations of the random network show that these values can strongly differ 
between realizations, and seem to be determined by the maximum value of di in a given 
network. 

As shown in fig. the phase <pi is negative for all distances and its absolute value in- 
creases with di. This implies that the response of oscillators to the perturbation is increasingly 
retarded as the distance grows. For small distances we find a linear regime, with a constant 
phase shift Aip between oscillators whose distances differ by one. Consequently, the pertur- 
bation propagates at constant speed v = \Cl/Acp\. For larger distances, before the regime of 
exponential decay in \Ai\ breaks down (cf. figs.^Ji and b), this linear dependence is lost. After 
an intermediate zone of faster variation, where the dispersion of the phase as a function of 
the distance is larger, the dependence of ipi on di becomes much less pronounced. The linear 
regime is particularly well defined for large perturbation frequencies, while the large-distance 
behaviour is clearly seen for small f2. 

The dependence of A4 on di can be understood by relating the connection matrix J with 
the distribution of distances in the network. The starting point is the solution to eq. (@J, 
A = Ai~ 1 e\. Since the eigenvalues of the connection matrix J are all less than or equal to z 
in modulus [6], the matrix Mr 1 = [(z + ifi)X — J]^ 1 can be expanded as a series of powers of 
J (for f2 ^ 0). The amplitude of the deviation from the fully synchronized state for the i-th 
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oscillator turns out to be 



where jf™' is an element of the m-th power of the connection matrix: J m 



A, = 



m=0 



(z + ift) 



-m-l t(«i) 



(5) 



{jg n) }. The 



matrix J m bears explicit information about the metric structure of the network. Specifically, 
the element jjj equals the number of directed paths of length m starting at node j and 
ending at node i [7]. Taking j = 1, i.e. the node at which the external perturbation is 
applied, we have in particular that = for m < di and j\™ ^ for m = di. If m > di, 

jjx^ is zero if no path of length m joins oscillators 1 and i, and positive otherwise. 

For oscillator i, thus, the first contribution to the sum in eq. |J5J is the term with m = di. 
If z/N <C 1 and di is small, there is essentially only one path of length di from 1 to i. This is 
realized by noticing that the probability of having more than one path of small length between 
two given oscillators is, at most, of order z/N. Consequently, for the majority of nodes at 
small distances from the perturbed oscillator, we have Jjf = 1. The total fraction n(di) of 
nodes with small di is also small, n(di) w z di /N and, therefore, the possibility that a node 
at distance di is also connected by a path of length slightly larger than di can be neglected. 
This implies that, for oscillators at small distances di, the sum in eq. (J5J is dominated by the 
first nonzero term: 



Ai ks (z + ift) 



-di-l 



(z 2 + n 2 y 



exp 



-i(di 



i ft 

l)tan _1 - 

z 



(0) 



This dominance is enhanced as z+ift becomes large, since higher-order terms are weighted 
by increasing powers of this number. In particular, the contribution of the first nonzero term 
becomes increasingly important as the frequency ft grows. In the right-hand side of eq. iJBJ, 
the exponential dependence of \Ai\ and the linear dependence of tpi on di, with a phase shift 
A(f = tan (ft/z), are apparent. Full straight lines in fig. ^show the quantitative agreement 
of this prediction with numerical results for small distances. Note that, as expected, the 
exponential approximation improves for larger ft. 

For large distances, such that z di ~ N, it is not any longer possible to insure that the first 
contribution to Ai is given by only one path of length di, nor that higher-order contributions 
are relatively negligible. In fact, one can argue that the number of paths of length m ending 



at a given node i scales as J 41 



(m) 



for large m. This can be proven inductively, noting 



that this number is z times the number of paths of length m — 1 ending at the oscillators to 
which i is coupled. This result is confirmed by the calculation of the average value of as 
an element of J m , under the hypothesis that the elements Jij of J are independent random 

variables. Assuming, thus, J-™' = Jqz" 1 , where the constant Jo may sensibly depend on the 
specific realization of the random network, we find 



J ) (z + ift) 



rn—di 



-m-l m _ I 1 i 



exp 



7T _i ft 

-l I — h di tan — 

2 z 



(7) 



Comparing with eq. JBJ, we realize that the dependence of \A+\ on di is now much slower, 
especially, for small ft. On the other hand, the phase shift between oscillators whose distances 
differ by one, Aip = tan _1 (ft/z), is the same as for small distances. When ft is small, the 
only difference is an additional shift of —n/2. This is clearly confirmed by numerical results: 
the dotted lines in fig. ^> are displaced by — n/2 with respect to the corresponding full lines, 
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Fig. 2 - Amplitude modulus at the maximal distance, di — 15, for the 10 3 -oscillator random network 
considered in fig. as a function of the perturbation frequency £1. The dotted line has slope — 1. 



which stand for the small-di analytical prediction. An independent confirmation of eq. J7J) is 
given by the fact that for small f2, large di, and a given connection network, the amplitude 
modulus should behave as ~ fl -1 . This is shown in fig. |21 for the same network as in 
fig. n at the maximal distance, di — 15. Note that this behaviour can be interpreted as 
a resonance-like phenomenon in the response of the oscillators -whose natural frequency is 
uji = 0- to the perturbation. The external harmonic action induces a larger effect when its 
frequency is closer to the individual frequency of the dynamical components of the system. 

Though the results discussed so far apply to random oscillator networks, our analytical 
approach can be used for arbitrary connection patterns. In the special case of regular networks 
we can even obtain exact results. Consider, for instance, a directed ring, i.e. a linear array 
with periodic boundary conditions, where each oscillator is coupled to its nearest neighbour 
to the left (z = 1). In this situation, Jy = 1 if i = (J + 1) mod N and otherwise, which 
implies that vanishes unless m = di + kN (k = 0, 1, 2, . . . ). Explicitly calculating the 

amplitude from eq. (JSJ, we find 

Ai = [1 - (1 + HI)-*] _1 (1 + ifi)"*- 1 . (8) 

Here, the amplitude modulus decays exponentially for all di and the phase shift between 
oscillators whose distances differ by one is constant, Aip = tan -1 n. 

Another situation where the decay of the amplitude modulus is purely exponential is the 
case of a tree-like connection network. Trees are graphs which contain no cycles. In directed 
trees, consequently, there is at most one path joining any two nodes, and therefore z = 1. 
Under such conditions, only one term in eq. JSJ) contributes to the amplitude, and 

A i = (i + m)- d *- 1 . (9) 

Those nodes which cannot be reached from the perturbed oscillator have Ai = 0. 

If the number z$ of connections ending at each node is not uniform over the network, eq. 
J2J can still be reduced to Q. Now, however, M = Z + ifll — J where Z is the diagonal 
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Fig. 3 - Amplitude modulus as a function of the distance from the perturbed oscillator in a 10 - 
oscillator random network where the number of neighbours of each node is Zj = 1, 2 or 3 with equal 
probability. Symbols are as in fig. The inset shows the amplitude modulus at the maximal distance, 
di — 15, as a function of the perturbation frequency. Axis scales in the inset are the same as in fig. 
[5] The dotted line has slope —1. 



matrix with elements Zij = In this situation, we cannot insure that the eigenvalues of 

the connection matrix J arc bounded in modulus, and therefore wc arc in general not able to 
expand M.~ l in order to obtain an expression similar to eq. JHJ. Whereas we may numerically 
solve eq. (0J by inverting Ai, our analytical approach does not apply anymore. Moreover, 
for such an arbitrary connection pattern the stability of the fully synchronized state is not 
guaranteed, and has to be separately analyzed for each realization of the network. 

Figure |3| shows numerical results for the amplitude modulus as a function of the distance 
from the perturbed oscillator, for a 10 3 -oscillator random network with a random distribution 
of Zi. Specifically, Zj is chosen to be 1, 2, or 3 with equal probability, in such a way that the 
average value of Zj equals the value of z for the network considered in figs. ^and^ Once Zj 
has been determined for each node, its Zj neighbours are selected at random. For the network 
realization of fig. the maximal distance is, again, di — 15. We find that, qualitatively, 
\Ai\ shows the same dependence on di and f2 as for the regular graph. As may have been 
expected, however, there is a larger dispersion in the amplitude modulus for elements at a 
given distance from the perturbed oscillator, especially, at large frequencies. The value of \A\ 
at the maximal distance, shown in the inset of fig. |31 exhibits the same dependence on the 
perturbation frequency as in the former case. 

In summary, we have shown that the response of a network of synchronized phase oscil- 
lators to a local harmonic perturbation of external origin exhibits two well-defined regimes, 
depending on the distance to the point at which the perturbation is applied. At small dis- 
tances, the system behaves as a dissipative medium. For a variety of connection patterns, 
the perturbation propagates through the network at constant velocity, while its amplitude 
decreases exponentially with the distance. In random networks, this behaviour holds up to 
certain threshold distance, beyond which the variation of the amplitude with the distance 
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is much slower. The smaller the perturbation frequency, the larger the amplitude at long 
distances and the smaller the threshold. Our analytical approach, which applies to regular 
graphs, suggest that the existence of these two regimes are a consequence of the geometrical 
properties of the network, associated with the number of paths connecting the perturbed os- 
cillator with any other node at a given distance. The analysis is linear, and the solutions for 
harmonic perturbations studied here can be combined to describe any other perturbation of 
sufficiently small amplitude. Nonlinear effects in the propagation of perturbations, which are 
expected to play a role for stronger perturbations, will be the subject of a future contribution. 

Let us finally point out that the present results establish a direct link between the dy- 
namics of an ensemble of coupled oscillators and the geometry of the connection network. 
This suggests a method to explore the structure of connections if the individual activity of 
oscillators is accessible through measurements. Determining the response of each oscillator to 
external perturbations applied to different elements of the ensemble could make it possible to 
reconstruct the underlying interaction pattern. 

* * * 

Valuable discussions with G. Abramson, H. Kori, M. Kuperman, and A. S. Mikhailov and 
are gratefully acknowledged. 
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